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INTRODUCTION 

It  is  proposed,  and  an  analysis  is  developed,  to  design  a  non- 
axial ly-symmetr ic  duct  intended  to  shroud  a  propeller.  One  of  the  first 
experimental  indications  of  the  possible  benefits  to  be  secured  from 
ducts  having  non-axially  symmetric  distributions  of  camber  was  given  by 
Shpakoff  and  Turbal  at  the  Twelfth  ITTC  in  Rome,  1 969 1 •  Unfortunately, 
the  official  proceedings  give  only  three  curves  which  compare  the  varia¬ 
tion  in  axial  blade  bending  moment  at  the  root  of  a  model  propeller  operating 
behind  a  model  of  a  tanker  with  a)  no  duct;  b)  symmetrical  duct;  c)  non- 
symmetrical  duct.  The  results  show  that  the  bending  moment  variations 
were  much  larger  with  the  symmetrical  duct  than  with  no  duct,  and  were 
greatly  reduced  by  the  duct  with  c i rcumferent ial 1 y  distributed  camber. 

This  experience  surely  refutes  the  oft-repeated  conception  that  a  duct 
homogenizes  the  inflow.  This  trend  is  confirmed  by  calculations  employ¬ 
ing  a  program  developed  by  Tsakonas,  et  al2  which  showed  that  vibratory 
thrust  was  predicted  to  be  increased  by  20  percent  by  a  simple  duct  as 
compared  to  the  ductless  case. 

In  1973.  Turbal3  advanced  a  theoretical  method  for  calculating  the 
shape  of  ducts  which,  in  his  terms,  "rectifies"  the  wake  flow.  Turbal 
attempted  to  calculate  a  duct  which  annuls  the  total  axial  and  tangential 
velocity  nonuniformities.  This  is  unnecessary  if  one  desires  merely  to 
reduce  the  blade  rate  shaft  thrust  and  torque  as  well  as  the  transverse 
and  vertical  forces  and  moments  about  the  horizontal  and  vertical  axes. 

It  is  also  unnecessary  if  one  merely  wishes  to  inhibit  the  onset  of  in¬ 
termittent  blade  cavitation.  With  the  first  objective  in  mind,  it  is 
only  necessary  to  calculate  additive  camber  distributions  which  induce 
opposing  spatial  flows  inside  the  duct  at  the  n-1,  n  and  n+1  harmonics 
(n  being  the  blade  number).  With  the  second  objective  in  mind,  it  is 
only  necessary  to  calculate  additive  camber  distributions  which  induce 
spatial  flows  opposing  the  lower  order  harmonics,  say  1,  2  and  3- 
Since  the  Fourier  harmonics  are  linearly  independent  and  the  theory  pre¬ 
sented  is  a  linear  one,  it  is  possible  to  design  a  duct  to  achieve  both 
objectives  simultaneously.  Turbal 's  analytical  method  does  not  allow 
any  insight  into  the  physics  of  the  problem  and  appears  to  suffer  from 
numerical  instabilities. 
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In  the  next  section,  an  analytical  method,  using  linearized  theory, 
is  presented  for  calculating  the  circumferential  camber  distribution  of 
a  non-axially  symmetric  duct  and  to  assess  the  amount  by  which  the  axial 
velocity  which  is  present  in  the  absence  of  the  duct  is  reduced  by  the 
duct . 

THEORY  OF  NON-AXIALLY  SYMMETRIC  DUCTS 

We  imagine  a  duct  without  a  propeller  immersed  in  the  wake  of  a 
single  screw  ship  or  submarine,  and  limit  our  attention  to  wakes  which 
have  port-starboard  symmetry  (i.e.,  the  ship  is  on  a  straight  course; 
the  submarine  could  be  trimmed  up  or  down,  but  not  yawed).  We  seek  to 
calculate  a  duct  camber  distribution  which  will  induce  axial  flow  compo¬ 
nents  opposing  those  of  the  given  wake  over  the  locus  of  the  disc  of 
the  propeller  at  spatial  frequency  m.  (The  frequency  m  may  take  on 
values  of  1,  2,  n-1,  n,  n+1 ,  etc.)  The  camber  distribution  must  be  such 
that  the  kinematic  boundary  condition  on  the  duct  surface  is  satisfied. 

The  unknowns  are  the  camber  distribution,  as  well  as  the  pressure  distribu¬ 
tion  on  the  duct,  as  functions  of  the  chordwise  and  circumferential  co¬ 
ord  i nates . 

We  begin  by  expressing  the  pressure  field  in  cylindrical  coordinates 
in  terms  of  a  pressure  jump  on  the  cylindrical  surface  s  =  r  "  : 


i  (x , r , y)  = 


2t 

/ 

o 


c 

/ 

-c 


fs  (R)sdSd9 


(I) 


where  R  is  the  Descarte  distance  between  the  fixed  point  (x,r,y)  and  the 
dummy  point  (£,s,@).  In  other  words 


R  {  (x-5)2  +  r2  + s2  -  2rs  cos  (9-y)  } ' ^ 


(2) 


By  expanding  1/R  into  a  Fourier  series  in  9-y,  it  can  be  shown  that 

2ir  c 
'm 

m-o 


P(x,r,y)=i  l  £m  /  /  ^17  [~  Vl/2((~"?)22rsr2  +  52)3cOS'n(Y'9)5d.Sf 

-  o  -c  /r  Vs  (3) 


Because  the  only  value  of  s  in  which  we  are  interested  is  s  =r^3  the 
two  symbols  will  be  used  interchangeably. 


where  e  =  I  for  m  =  0 
m 

=  2  otherwise 


Q  .  is  the  associated  Legendre  function  of  the  second  kind  of  half 
m-  1  /  / 

odd  integer  order.  The  properties  of  the  Legendre  function  are  given, 
for  example,  in  Ref.  4. 

Upon  carrying  out  the  indicated  differentiation  with  respect  to  s, 
the  field  pressure  becomes 


i  e  271 

I  r  m 


P=^tI—  /  cosm(y-e)d8  /  Ap  (e  ,  C)  I  md^ 


>'rs  o 


where 


4rs  (m  +j)  [s2  -  r2  -  (x  -  £)2] 

=  - *- 


z  0 


1/2  [(x-02  +  (r-s)2][(x-5)2+  (r+s)2]  m+,/2  m"1/2 


(4) 


(5) 


and  the  argument  of  all  the  Q  functions  is 
(x-S)2  +  r2  +  s2 


Z  = 


2rs 


(6) 


If  the  pressure  jump  Ap  is  expressed  as  a  Fourier  series  in  8, 


-P  =  l  iPn(?)cosne  (7) 

n=o 

and  the  field  pressure  is  also  expressed  as  a  Fourier  series  in  y. 


X> 

p  =  y  P  (x , r) cosmy 
m 

m=o 


(8) 


Then,  upon  substituting  (7)  into  (4)  and  carrying  out  the  8-integral, 
it  can  be  shown,  as  a  consequence  of  the  orthogonality  between  the  trig¬ 
onometric  functions,  that 


When  r  =  s  ana  •,  =  x,  Z  becomes  unify  2nd  the  functions  are  to-jarithmiczl 
infinite.  Although  this  singularity  is  integrable ,  it  offers  problems  in 
the  numerical  evaluation  of  the  integrals  which  can  be  avoided  altogether 
taking  the  maximum  value  of  r  to  be  s-e  where  c  z>  a  small  quantity. 
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Pm(x*r)  = 


e  c 
m 


/ ap  (e)  id? 

_  m  m 

8tt  t'TT  c 


(9) 


As  a  consequence  of  the  linearized  Bernoulli  equation 


3<j> 

P  =  PU^ 


(10) 


Equation  (9)  can  be  restated: 


,J_  _  m 

U  3x 

m  IbTiVrs  -c 


/  ACP  (?)  I  <U 

'  m  m 


(11) 


where  is  a  pressure  coefficient.  Equation  (11)  gives  a  relationship 
between  the  axial  flow  inside  the  duct  and  the  pressure  jump  across  the 
duct  sections  for  every  Fourier  component.  If  the  pressure  jump  were 
known,  the  axial  velocity  could  be  calculated.  It  is  our  desire  to  choose 
a  pressure  jump  such  that  the  induced  axial  velocity  nullifies,  as  much 
as  possible,  the  nonuniformity  of  one  or  more  Fourier  components  of  the 
ship  wake  and  that  has  shockless  entry  at  the  leading  edge.  Thus,  Eq. 

(11)  is,  in  fact,  a  Fredholm  integral  equation  of  the  first  kind  to  de¬ 
termine  X  in  terms  of  (l/U  '-i/'-x)  .  Examination  of  the  corresponding 
Pm  m 

two-dimensional  integral  equation  has  disclosed  that  the  solution  is  not 
unique.  Presumably,  the  three-dimensional  equation  given  here  has  the 
same  property.  If  so,  the  property  is  advantageous  because  it  means 
that  there  is  more  than  one  pressure  distribution  that  will  do  the  job. 

The  following  method  was  originally  proposed  to  solve  Eq.  (11): 

We  define  a  rooftop  pressure  distribution  for  the  harmonic  as 

fo 1 1 ows : 


AC  (■:)  =  h, 
P 


-  r  < 


x  (-:)  =  h (— — — ),  x 

p  c  -  X  ’  o 

o 


(12) 


The  lift  coefficient  on  the  section  is  then 
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C 


i 


2C 


/ 

-c 


AC  d£ 

P 


[3  +  ^3 


(13) 


When  h  =  I,  the  resulting  distribution  is  said 
now  consider  several  unit  rooftops  that  differ 
of  xq  and  distinguish  them  using  the  subscript 
distribution  will  consist  of  a  linear  sum  of  M 
weighted  by  unknown  coefficients  a..  Thus, 


to  be  a  unit  rooftop.  We 
only  in  the  numerical  value 
i .  The  correct  pressure 
of  these  unit  rooftops 


AC 

P 


M 

y  a .  AC 

.1  '  P 


(14) 


Corresponding  to  each  unit  rooftop  is  a  radial  distribution  of  axia  city 

at  a  specified  station  x  designated  as  u..  Thus 

u  (r)  =  T a.u.  (r)  (15) 

c  i  i 


where  the  subscript  c  means  that  the  velocity  is  calculated.  We  now  con 
sider  a  penalty  function 


I 


s-e  2 

[  rdr  [u  -  u  ] 
i  me  c 

o 


(16) 


where  the  subscript  me  means  that  the  velocity  is  measured.  The  weighting 

factor  r  is  introduced  to  emphasize  the  outer  radii  of  the  duct  where  the 

propeller  is  most  heavily  loaded.  We  wish  to  determine  u^  so  that  it 

equals  u  as  best  it  can.  We,  therefore,  minimize  I  with  respect  to  each 
m 

of  the  a.'s  and  set  the  minimum  equal  to  zero.  Upon  substituting  (15)  into 
(16)  and  carrying  out  these  operations,  we  are  led  to  the  following  system 
of  equations  for  the  a.'s. 


s-e 
/  rdr 
o 


u  .  u 
J  m 


M  s-e 

y  a .  j  u.u,  rdr 

,  i  ii 


I  ‘  j  <  M 


(17) 


s-c 

The  coefficients  J  u^u.rdr  forma  symmetric  matrix.  Calculations  of  the  u.'s 
and  the  mat  r  i  x  coe?  ficients  werecarriedout  using  Eq.  (11)  for  five  different 
unit  roo  ft  ops  with  x=0andx  =-1.0,  -0.6,  -0.2,  +0.2,  +  0.6  leading  to  a  five- 
by-five  matrix.  It  was  found  that  the  matrix  was  highly  ill-conditioned. 

This  means  that  the  axial  flow  in  the  duct  is  unable  to  distinguish  one 
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unit  rooftop  from  another  or,  put  another  way,  that  the  shape  of  the 
axial  velocity  pro^le  in  the  duct  is  virtually  the  same  regardless  of 
the  shape  of  the  pressure  distribution  on  the  duct.  For  this  reason, 
the  concept  of  representing  the  pressure  distribution  on  the  duct  by 
several  unit  rooftops  had  to  be  abandoned.  It  was  recognized  that  the 
velocity  profile  in  the  duct  is  primarily  dependent  on  the  magnitude  of 
the  total  sectional  lift  and  is  insensitive  to  the  shape  of  the  distribu¬ 
tion.  Thus,  we  selected  only  one  unit  rooftop  to  represent  .iC  •  The  one 

P 

chosen  was  for  =  +0.6  and,  in  this  case,  the  matrix  equation  (17)  re¬ 

duces  to  a  single  equation  for  a.  = h.  The  result  is 

s-e 

j  r dr  u  u, 

J  me  * 

h  =  -2— -  (18) 

rdr  u,- 
o 

Once  h  is  determined  in  this  way,  the  lift  coefficient  on  the  duct  section 
can  be  calculated  using  Eq.  (13). 

It  is  stated  on  page  529  of  the  paper  by  Van  Manen  and  Oosterveld3 
that  the  sectional  lift  on  an  axially  symmetric  duct  should  not  exceed 
unity  at  the  risk  of  inducing  stall.  In  our  case,  fhe  lift  varies  circum¬ 
ferentially,  alternating  between  positive  and  negative  values  since  it  is 
proportional  to  cos m  . .  Because  the  maximum  lift  is  only  experienced 
local Iv  as  we  proceed  around  the  c  i  rcu-fererce .  the  Van  Manen  and 
Oosterveld  limitation  may  be  too  restrictive.  However,  it  is  certainly 
a  conservative  limitation  and,  in  the  absence  of  any  other  criterion,  we 
will  use  it. 

Mow,  the  lift  coefficient  calculated  using  Lqs.  1181  and  (13)  may 
be  smaller  or  greater  than  unity.  If  it  is  greater  than  unity,  it  must 
be  reduced  to  unity  in  accordance  with  our  criterion  and  the  value  of  h 
must  be  correspondingly  reduced.  This  means  that  the  duct  can  be  expected 
to  cancel  the  m*^  harmonic  of  the  measured  axial  velocity  only  partially. 

We  present  now  an  example  of  how  a  duct  may  be  expected  to  affect 
the  axial  flow.  For  the  example,  only  the  first  harmonic  of  a  wake  will 
be  considered.  The  data  is  taken  from  Figure  27  of  hadler  and  Chengf 
which  is  reproduced  here  as  Figure  1.  Only  the  data  for  the  first  harmonic 
and  Model  4602,  indicated  by  circles,  will  be  considered.  The  chord  of 
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the  duct  is  taken  to  be  equal  to  the  radius  (s=2);  the  section  where 
the  first  harmonic  is  to  be  rectified  is  the  midsection  (x = 0) ;  the  roof¬ 
top  distribution  begins  to  fall  off  at  0.8  of  the  chord  (x  =.6).  Then 

C  ,  as  calculated  using  Eqs.  08)  and  (13)  becomes  C„  =  -9766,  which  is 
l  1 

within  the  Van  Manen  and  Oosterveld  limitation.  A  plot  of  the  data  from 

Figure  1  is  reproduced  in  Figure  2,  together  with  the  modified  wake  harmonic 
to  be  expected  with  the  duct  in  place.  It  can  be  seen  that  the  large 
velocities  at  the  outer  radii  have  been  significantly  reduced.  At  the 
inner  radii,  the  negative  velocities  have  become  more  negative.  On  the 
other  hand,  it  is  well  known  that  intermittent  cavitation  occurs  primarily 
at  the  outer  radii  where  the  bulk  of  the  blade  loading  takes  place,  and 
so  the  beneficial  effects  of  reducing  the  axial  velocity  at  the  outer 
radii  should  more  than  compensate  for  the  deleterious  effects  of  increas¬ 
ing  the  axial  velocity  at  the  inner  radii. 

It  is  the  ultimate  purpose  of  this  project  to  design  and  build  a 
duct  according  to  the  precepts  set  forth  here.  Consequently,  a  series 
of  wind  tunnel  measurements  taken  behind  a  wake  screen  were  carried  out 
at  the  Princeton  wind  tunnel  by  Professor  R.  Hires.  The  screen  was  designed 
to  simulate  the  flow  in  a  ship  wake.  Details  of  the  experiments  will  be 
presented  in  a  subsequent  report.  For  present  purposes,  we  present  only 
the  result  of  the  amplitudes  of  the  first  harmonic  of  the  measured  wake  . 
This  is  shown  in  Figure  3.  The  results  of  carrying  out  the  computations 
indicated  by  Egs.  (18)  and  (13)  show  that  C  =3-0.  Since  this  exceeds 

V 

the  Van  Manen  and  Oosterveld  limit,  C  was  reduced  to  unity  and  the  modified 
wake  with  the  duct  in  place  is  also  shown  in  Figure  3-  It  can  be  seen 
that  the  modification  created  by  the  duct  is  a  smaller  percent  of  the 
total  in  this  case  than  it  was  for  the  Hadler  and  Cheng  wake  previously 
analyzed.  This  is  because  the  lift  coefficient  had  to  be  reduced  to 
unity  in  accordance  with  the  Van  Manen  and  Oosterveld  limitation. 

THF  GEOMETRY  OF  THE  DUCT  SECTIONS 

Up  to  this  point  we  have  examined  the  effect  of  a  duct  rooftop 
load  distribution  on  the  flow  inside  the  duct.  At  this  juncture,  we 

*  The  rhase  varies  slightly  from  radius  to  ra-iius  even  though  the  wake 
screen  is  rresunablj  synrr.etric.  The  iuct  cannot  be  expected  to  follow 
this  radial  vzri  :t ion  in  phase. 
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Now  consider  the  linearized  boundary  condition  on  the  duct.  This 
'  <-  that  I  /l-  •  ■:  /  -  r  *  slope  of  the  duct  at  r  =  s.  Now,  the  slope  consists 

C  two  parts,  one  due  to  the  pitch,  >  ,  of  the  section,  and  the  other  due 

rn 

: o  t u  e  slope  of  the  camber  line.  Thus,  if  the  camber  line  is  defined  as 

\  \  ,  :  ”e  fixmcai  \  condition  is 

'r  -  i  +  V Jx)  (25) 

_hpn  Eg.  T  3  •  becomes 

v  '  ,  x  f  ■=  -i  o  i,  s )  (26) 

Upon  integrating  with  resnect  to  \  from  the  leading  edge  where,  bv  defini¬ 
tion,  vtp)  =  0,  we  obtain 


v  t  s  ^  =  .  t  x  +■ 1  '  *■  g  t  s  '  )  d  \ 

in  i"i  •  "i 


(27) 


We  define  v  to  be  zero  at  the  trailing  edge  also,  so  that,  upon  setting 

ni 

<  -  *1,  we  obtain  the  following  expression  fo*  , 

tTt 


2 


I 

i 

-1 


q  (  x  '  )  d  x  1 

■it 


(28) 


Substituting  for  g  1  *  '  )  from  f.q.  (23)  and  interchanging  the  order  of 
integral  ion 


/?  I 

-ti 

“ 


k  /  ■  ).|- 


(29) 


(30) 


‘-;e  -|e f  i i  t  i  on  >‘  Eg.  (24),  and  letting  u  *  t  “  X,  we  find 


}y  H  ( u )  d  u 


(31  ) 
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This  double  integral  can  be  manipulated  into  two  single  integrals  as 
shown  in  Appendix  A.  The  result  is 


K  U)  =  /  (1+u-C)  H  (u)du  +  2  /  H  (u)du  (32 

m  -(1-0  m  1+S  m 

The  point  u=0  must  be  excluded  because  the  integrand  is  singular  in  the 
Mangier  sense.  The  final  result,  ready  for  numerical  computation,  is 


K  (c)  =  f  ( 1  't'u* OH  (u)du+  /  (l+u-t;)H  (u)du  +  2  /  H  (u)du - (1-0 


-0-0 


where  c  is  a  small  quantity. 

To  summarize  the  computation  of  a  ,  it  is  determined  from  Eq.  (29) 

m 

with  K ( O  given  by  Eq.  (33)  and  H^(u)  given  by  Eq.  (20)  with  Z  =  l+u2/2s2. 
Moreover,  the  kernel  K(0  is  logarithmically  singular  at  £  =  ±  1  and,  al¬ 
though  these  singularities  are  integrable,  they  present  problems  in  numer¬ 
ically  calculating  the  integral  in  Eq.  (29).  The  problems  can  be  avoided 
by  changing  the  limits  on  the  integral  in  Eq.  (29)  from  ±1  to  ±(l-e)  where 
c  i s  a  sma 1 1  quant i ty . 

In  a  similar  manner,  from  Eq.  (27)  we  can  obtain  an  expression  for 

the  ordinates  of  the  camber  line  with  a  considered  to  be  known. 

m 

£  /2  1 

y  -  a  (1+x)  =  -7T—  /  AC  (0d€  K  (0  ( 3 k) 

m  m  rtTT  ' .  n  x 


where 


Kx  (0  =  /  Jm(xrcjdxi 
m  - 1 


or,  using  the  definition  of  J  ,  Eq.  (24) ,  and  letting  u  =  t-x 

m 


Kx  (O  =  /  dx  /  H  (u) 
m  -I  '^-xl 


In  a  manner  similar  to  that  presented  in  Appendix  A,  the  order  of  integra¬ 
tion  can  be  interchanged,  the  x^  integral  can  be  carried  out,  and  the 
final  result  is 
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K  =  /  H  (u) [x+u~c]du  +  (1+x)  /  H  (u)du 

xm  5-x  m  1+C  m 

-e  1+C 

K  =  f  H  (u)[x+u~C]du  +  /  H  (u)  [x+u~C]du 

x  t  m  '  m 

m  -(x-C)  e 


C  >  x 


+  (1+x)  /  H  (u)du  —  (x-C) 

1+C  m  e 


C  <  x 


where,  again,  e  i s  a  small  quantity. 

To  summarize,  the  ordinates  of  the  camber  line  are  given  by  Eq.  (3*0 

where  Kx  is  defined  by  Eqs.  (37)  and  (38)  and  H  (u) ,  Z,  a  are  defined 
m  mm 

as  before.  Once  again,  the  kernel  is  logarithmically  singular  at 
C  =  -l  and  C  =  x.  The  singularity  can  be  avoided  by  changing  the  lower 
limit  in  Eq.  (34)  to  -(1-e)  and  surrounding  the  point  C  =  x  by  a  small 
interval  which  is  excluded  from  the  integration  scheme. 

THE  RESULTS  OF  SOME  NUMERICAL  COMPUTATIONS * 

The  value  of  a  for  several  values  of  m  were  calculated  for  the 

m 

case  s  =  1  ,  x  =  0,  xq  =  0.6.  The  results,  together  with  one  isolated  case 
with  s=2,  are  shown  in  Table  1. 


Table  I.  Values  of  a 

m 

a  /C  (s-l)  a/C  (s-2) 

mi.  m  / 

1  0.2912  0.1563 

2  0.5607  N. A. 

3  0.8322  N. A. 

4  1.1038  N. A. 

5  1.3356  N.A. 

It  can  be  seen  that,  as  m  increases,  the  angle  a  increases,  and 

m 

this  means  that  the  duct  section  must  work  harder  and  harder  to  produce 
the  necessary  axial  flow  to  cancel  a  unit  wake  inflow.  Fortunately,  the 
magnitude  of  the  wake  inflow  as  the  order  of  the  harmonic  increases  tends 

A  subroutine  for  calculating  the  Legendre  functions  was  developed  based 
on  a  series  given  in  Reference  4  for  1  <  Z  <1.09  and  based  on  a  standard 
hyper geometric  series  for  Z  2  1.  09  , 
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to  get  smaller  and  smaller  so  that,  in  practice,  the  higher  harmonics 
can  usually  be  dealt  with  without  exceeding  the  Van  Manen  and  Oosterveld 
1 im i t . 

The  cases  m=l,  s  =  1  ,2  indicate  that  a  greater  value  of  C?  is  re¬ 
quired  for  the  case  s=2  (chord  -  rad i us)  than  for  the  case  s  =  1  (chord  = 
diameter)  to  achieve  the  same  axial  flow  in  the  duct,  and  this  is  intuitively 
correct. 

Points  on  the  camber  line  for  the  case  m=l,  s=2  were  calculated 
and  are  shown  in  graphical  form  in  Figure  A  and  in  tabular  form  in  Table 
II.  These  points  were  fitted  to  a  polynominal  passing  through  the  points 
(1,0)  and  (-1,0)  using  a  least  square  fit,  and  the  resulting  polynomial  is 


~  =  O'*2)  [a  +  a  x  +  a  x2  +  a  x3  +  a  x4] 
0  0  12  3  4 


where  a  =0.15577 
o 

aj  =  -0.01669 

a2  =  0.05212 
a3  =  -0.055A0 
a4  =  -0.03812 


(39) 


It  should  be  pointed  out  that,  because  we  are  dealing  here  with  the 
first  harmonic, Eq.  (39)  must  be  multiplied  by  cos0 ,  and  also  that  a,  must 
be  multiplied  by  cos0.  Thus,  both  the  camber  and  pitch  angle  vary  circum¬ 
ferentially.  Since  the  sections  are  at  ijcose,  this  effect  can  be  achieved 
simply  by  tilting  the  complete  duct  at  angle  .  For  all  higher  harmonics 
no  such  simple  recourse  is  available,  and  the  local  pitch  angle  must  be 
built  into  the  duct . 


It  is  because  of  the  relative  simplicity  of  the  first  harmonic  that 
the  Princeton  tests  and  the  computations  carried  out  thus  far  have  con¬ 
centrated  on  this  case.  The  ultimate  purpose  will  be  to  build  a  duct 
based  on  the  camber  line  of  Eq.  (9),  to  tilt  it  at  the  angle  otj  =  .  1 563 , 
according  to  Table  I,  and  to  mount  it  in  the  Princeton  tunnel  behind  the 
same  wake  screen  as  used  in  the  preliminary  tests  and  then,  by  taking 
measurements  inside  the  duct,  to  prove  the  efficacy  of  a  non-axially 
symmetric  duct  by  demonstrating  a  reduction  in  the  first  harmonic  of  the 
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wake  as  shown  theoretically  in  Figure  3-  Such  a  duct  is  shown  sche- 
matical 1y  in  Figure  5- 

CONCLUSIONS 

A  theory  for  the  behavior  and  design  of  non-axially  symmetric  ducts 
intended  to  shroud  a  propeller  is  presented.  Measured  wake  data  is  also 
presented,  together  with  the  modification  of  the  wake  caused  by  an  appro¬ 
priately  designed  duct.  Moreover,  the  shape  of  such  a  duct  is  calculated 
for  one  case.  It  remains  to  build  a  model  of  such  a  duct  and  measure 
the  modified  wake  inside  it  to  demonstrate  the  tenets  of  the  theory,  and 
also  to  show  that  non-axially  symmetric  ducts  can  be  expected  to  be  useful 
in  ameliorating  transient  propeller  cavitation  and  vibratory  shaft  forces. 
The  results  of  this  ongoing  program  will  be  presented  in  a  subsequent 
report . 
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TABLE  II* 

POINTS  ON  THE  CAMBER  LINE  OF  THE  DUCT 
TO  BE  USED  IN  PRINCETON  WIND  TUNNEL  TESTS 


X 

V/c, 

.9 

.0280361 

.8 

.0564023 

.7 

.0847943 

.6 

.109061 

.5 

.126248 

.4 

.138521 

•  3 

.146653 

.2 

.  153209 

0 

.156918 

-.2 

. 150040 

-.3 

. 143592 

-.4 

.134574 

-.6 

.107941 

-.8 

.0682754 

U?  .1  . 
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APPENDIX  A 

REDUCTION  OF  THE  DOUBLE  INTEGRAL 


1  °° 

K  (£)  =  /  dx  /  H  ( u ) d u 
m  ,  '  m 

-1  £-x 

In  the  u,x  plane,  the  integral  is  carried  out  over  the  shaded  re 
gion  defined  by  the  sketch  below: 


If  the  order  of  integration  is  changed,  we  must  integrate  on  x 
first.  Thus,  we  integrate  over  the  triangle  defined  by  the  lines  u  =  £  -  x 
u  *  5  +  1  ,  x=+l,  and  then  add  the  integral  over  the  sem  i  -  i  nf  i  n  i  te  rectangl 
def  i  ned  by  the  lines  u  =  ►  +  1  ,  u=^°,  x  =  :  I.  The  result  is 


(1  +  0 

K  U)  -  f  H  (u)du 
m  .  ■ ‘  .  m 


1  »  +1 

f  dx  +  /  H  (u)  f  dx 

:>u  +£ 


After  the  x-integrals  have  been  carried  out,  the  result  is  Eg.  (32). 


APPENDIX  5 
COMPUTER  PROGRAMS 

The  five  computer  programs  that  were  used  in  the  numerical  computa¬ 
tions  are  presented  herein,  together  with  descriptions  of  what  they  do 
and  with  a  dictionary  for  each,  relating  Fortran  variables  to  analytic 
variables  for  input  and  output  statements. 
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PROGRAM  1 

Calculates  a/C^  for  given 

FORTRAN  ANALYTIC 

VARIABLE  VARIABLE 

S  s 

M  m 

XO  x 

o 

A 

Input  B 

AERR 
RERR 

EPSK  e 

AERRK 
RERRK 


ALPHA  .i/C 

Output  CK 


DESCRIPTION 


Lower  limit  of  outer  integral 
Upper  limit  or  outer  integral 
Absolute  error  tolerance  for  outer  intrgral 
Re  1  at  ive  error  tolerance  for  outer  integral 
Lower  limit  of  certain  inner  integrals 
Absolute  error  tolerance  for  inner  integral 
Relative  error  tolerance  for  inner  integral 


Value  of  outer  integral 


« 


> 


B-2 
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EXTERNAL  ?  S  2 
C CM MOV  S,M,3 

>•  COMMON / COMO 1  /CO  C  5 3 , 3  )  ,  AQ  ( 3  )  ,  A2  Q  ,  3Q  ( 3  ) ,  GO  (  3  ) , 

1  A1Q(3),A3j(3),E0(3),S1U(3),33Q(3),FAC(3) 

CCKMC,V/COVQ2/ A222(6,3),BQ22(6,3  ) 

)  CCMMO.N/COVFT/DQ  (  3  ) ,  SO,  AM3 

C0MM3,!/C0'«,K2/X3,C,CL,E?SK,  AERRK,RER?K 
COM‘:ON/COMl/LCl,LC2,LC3,LC4,LC5 
LC1  =  D 
LC2  =  . 

LC3  =  C 
LC4  =  J 
LC5  =  J 

RZAD(5,9D1 )  3,M,XJ 
m  ?ORM  AT(F,I,F) 

RFAD(5,9J2)  A, 6, AE  RR, RE  RR 
9  32  F  ORV  A  7  (  4F  ) 

READ (5,932)  EPS F, ASRRK, RERRK 
CALL  ’JSRSET  ( l ,  LEVL ) 

SC=S*S 

r  1  =  3. 1415927 

• »  —  i  • 

C  =  l. 

CL=  1 . 

CALL  QCONST 

CK  =  FCADRE( FX2, A,  5, AERR, RERR, ERR , I ER  ) 
PX=1./(16.*PI*Q> 

AL?HA=PX*CK 
WRITE(6/9M3)  ALPHA, CK 

933  FCRMATdX,  'ALPHA  =  ',  E14 .6,  2X,  '  CX  =  ',514.6) 

END 


J 
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r  UNCTION  F  X  2  (  G  ) 

EXTERNAL  FT,FT1,FT2 

J  ccx.'cn  s,-,gi 

CCMXC?:/CjV(jl/C^(5?^/3)/AQC3)/A2Q/93C3)/GS(3)/ 
1  AlG(3),A3Q(3)/E3(3)/SlQ(3)/33u(3)/rAC(3) 

3  CuMMO‘J/COM.32/A322(6/3)/3Q22(6/3  ) 

CCMM0K/CCMFT/D9(3 )  , SG / A V3 
C CMM C  .*!/  COX  X 2  /  ;<  2  ,  C ,  CL ,  E?  S K  ,  A  E R RK  /  RE RR K 
CCXVCN/C0M1/LC1,LC2,  LC3/LC4/LC5 
LC5=L75+1 
G1  =  G 

ai  =  eps:< 

31=1 .-G 
A  2  =  EPS K 
52=1 . +G 
A 3  =  1 .  -t-G 
o3=12.5*S 

aerr=aerp:< 

RERR=?ERRK 
F  e  1  =  D  . 

IF  (C.EG.l .  )  G3  :C  11 

Cl  =  rCArFE(r71/Al,cl/AERR,EER  R/ERROR/IES) 
C2=DCADFZ(FT2/A2/S2/AERR/RERR/£RROR/ISR) 

C3=DC ADPX( r 7/A3/S3/AERR/RERR /ERROR/ IER) 
F71=C1«-C2>2.*C3-4.*(1.-C)/EPSK 
11  CONTINUE 

H=4.*CL/(3.+X3/C) 

IF(-C. LT.G. AND.G.LS.X3)  DP  =  H 
IF(XB.L7.G.  AND.G.L7.C)  DP=H*  ( C- G  )/  ( C-XU  ) 

F  K2  =  F  K 1 *DP 

RETURN 

END 


o 


D 


B-b 

* 
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e  U  N  C  T  I  a  N  T  (  T  ) 

COMMON  S/M/GG3G 

J  C  C  M  V  0  N  /  C  0  M  ;  ?  /  D  i(3),  5  .rA‘-V 

COM  MO  ?;/  CO  !  1 / LC1 ,LC2 , LC3 , L  C4, LCd 
LC2  =  LC2  +  1 

~)  DEL=T*T/(2. *J3) 

Z=DE0+1. 

Z1=Z*1. 

ZQ=Z*Z 

DEL2=DEL*DEL 

IF(DZL.GE.H.jy)  GD  TO  11 
IF(DZL.LT. l.E-7)  GO  TO  9 
G  1  =  011  2  (  DEL/  1  ) 

U2=GG2(P£L, 2 ) 

23  =  0 , 2 ( DEL, 3  ) 

iuiP.-il  =  D  3  (  2  )  /  (2.+0EL  )  *  (  0  2  -  (  1 .  *  D  Z  L  )  *  Q 1  ) 

GPk2  =  Du(J)/(2.OEL)*(0  3-  (l.+L>EL)*G2) 

FT  =  l./(SG*(2.+DEL))*((2.«,DEL)/2.*Ol  + 

1  2.*DG(2)*((-2.+L'rL)*G2<-(2.+DEL-DEL2*(3.-*‘DEL))*Ql)/ 

1  (  (2.+-0EL)  *DEL)»QPK2-(  1  .  +  2EL)*jr>Kl  ) 

RETURN 

9  continue 

FT  =  1  . /(2.*33)*(  (r  V*F’1-.75)*ALCJ(DZL)*2./D-L) 

RETURN 

11  CON  T  I  NUE 

Gl  =  3Gl  ( 2, 1  ) 

02=301(2,2 ) 

Q 3  =  Q0  1  ( Z, 3  ) 

CPlsDQ(2)/Zl*(n2-2*Ul ) 

QP2=CQ(3)/Z1*( 33-Z*Q2 ) 

FT=l./(SQ*Zl)*(Zl/2.*Ql+2.*DO(2)*(((Z-3.)*u2- 
1  (ZQ*Z-4.*Z+1. ) *0 1 )/(Z3-l. )+QP2-Z*3?l) ) 

;  RETURN 

tl  *.f  n 
CC 

J  FUNCTION  F  T 1 ( T ) 

COMMON  S  ,  M. ,  G 

COMMON /CO’1  l/LCl,LC2,LC3,LC4,LCb 
O  LC1=LC1+1 

FT1=( l.-T-G) *FT( T) 

RETURN 

O  end 

CC 

FUNCTION  ?T2(T) 

COMMON  S ,  M ,  G 

CCMM0N/CDul/LCi,LC2,LC3/LC4/LC5 

LC1=LC1>1 

rT2=(l.+T-G)*ET(T) 

RETURN 
E  N  C 
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PROGRAM  2 


for  given  m,  3.  xq,  a/C^ 


ANALYTIC 

VARIABLE 


DESCRIPTION 


Cal cu I  ate  Y (x) /C 

FORTRAN 

VARIABLE 

s 

M 

XX 

XO 

Input  ALPA 

A 
B 

AERR 

RERR 

EPS 


s 

rn 

x 

x 

o 

a/C 

L 

Lower  limit  of 

Uppe r  limit  of 

Absolute  error 

Relative  error 

e  Integrate  from 

XX  +  EPS  to  B 


first  outer  integral 
second  outer  integral 
tolerance  for  outer  integral 
tolerance  for  outer  integral 
A  to  XX-EPS  and  from 


Output  YY 

ALPA 
CK 
CKI 
CK2 


Sum  of  1st  and  2nd  integrals 
1st  Integral 
2nd  integral 


l 


i 
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LXTErNAL  F  K 2 
COMMO'.’  S,A,Z 

CC-!'-,CN/Cu“-l  /  C  - j  C  t  .3  3,3),  AU(3  ),  A2C,3G(3  ),  (  3  3, 

1  A1C(3),A3C(3),E3(3),F10(3),E30(3),?AC(3) 

COMMON/ COM C.2/ A* 22  (6, 3),  SC 22  (5, 3  ) 
CC.VMCN/CO.MFT/D^U  ),S0,  AM? 

cc;'M.c?.7ca*-*x2/X3,c,CL 

COMMCV/COMl/LCl,LC2,LC3,LC4,LC5 
COM MON /COM XX/ XX 
LC1  =  3 
LC2  =  ? 

LC3  =  ? 

LC4  =  0 
LC5=2 

READ(5,9<)1)  3,  M,  XX/  XO,  ALPHA 
9t?l  r OKMAT(F/I,3f  ) 

READ (5,9/2)  A, 3, AER R/ KERR, EPS 
9/2  FORM A  T ( 5r ) 

CALL  UERS"?( 1,LE7L) 

SQ=S*S 

PI =3 . 1415927 

0=1 . 

A 1  =  A 

B1=XX-L?J 
A2=XX+EPS 
82=E 
C=1 . 

CL=1. 

CALL  QCONST 

CKl=r CADRE (FK2, A1 ,81 , A ERR, R ERR/ ERR, I ER ) 

WRITE  (  6/  9*3  J  33  )  CXI 
9311333  F0k.MA?(lX,El4.o) 

CK2=F CADRE (FK2, A2,E2, AERR,RERR, ERR, IER) 

ck=c;<i+ck2 

PX= 1 # / ( 8 . *  P I *0 ) 

YY=PX*CK-ALPhA*(XX+l. ) 

WRITE (5, 933)  YY, ALPHA, CK, CK1, CK 2 
9/3  F0RMATC1X,  '  Y  =  ' , El  4 . 6, 2 X ,  '  ALP n A  =  *,E14.5,2X, 

1  /IX, 'CK  =  ',3E14. o  ) 


r-2298 


S-  VICTIM*.’  TV  2(G) 

LXT-r.'IAL  r  T  ,  c  »T  1  ,  r  T  2 
j  C C V. l-i 2 ?!  3/  u  ,  G 1 

CO*4MO\’/C3MGl/C*(5G2,3)/AG(3)/A2Q,BQ(3),GG(3), 
1  AlG(3),A3v«(3),E'J(3),9lQ(3),33'H3),FAC(3) 

~)  CCMyGS’/CiJVQ2/A322(6,3  ),RQ22(6,3  ) 

COMVON'/COMr  T/ DG  (  3  ),  SO,  A ".3 
CCV.ya\/C3MK2/:<3/C,CL/EPSK 
Cny.MQNVCnMl/LCl/LC2/LC3,LC4/LC5 
ccv.vn:j/CQyxx/xx 

LC5=LC5+1 

1  —  o 

SI.\'rIN  =  12.5*3 
A  1 1  =  E  P  S  K 
tll=XX-G 
A 1 2  =  £  ?  S  K 
312  =  1  .  +•  G 
A  1 3  =  1 . ♦G 

e  1 3 = s  i  r:  f  i  \ 

A22=:-XX 

b  - 1  •  *  u 

A  2  3  =  1  ,  +  c 
a  2  3  =  S  I M  F  I 
AEHf sS. 

RE.RRs.tfWJl 
F  X  X  =  1 .♦XX 
FX1=3. 

IF  (G.EQ.l.)  GO  TO  11 
IF(G.GT.XX)  GO  TO  111 

C1  =  CCADPE(  FT1,  AU,B1 1,  AERR/RERR/ERROR,  IER) 

C2  =  DC  ACRE  (FT  2,  *12,t>12,  \ERR,RERR,ERROR,IEk) 
C3=DCAPRE(rT,A13,El3,AESR,RERR,ERP3R,IER) 

F F l=Cl*C2*r XX*C3-4. *(XX-G)/E?SK 
GO  TO  11 
111  CCMINUS 

J  Cl  =  CJA0OE(r'T2,  A22/022,  AERR, RERR , ERROR,  I  ER  ) 

C2  =  0CADRZ(?T,A23,323/ AERR, RERR, ERROR,  I ER) 
r K1=C1*FXX*C2 
O  11  CONTINUE 

»:*•». *CL/(3.»X3/C) 

IF(-C.LT.G.  AND  .  G.LE.  X«3  )  DP  =  H 
O  IFCX.’.LT.G.  AND.G.LT.C)  DR=H*(  C-G  )/ (  C-XG  ) 

r  K2  =  Fl'l  *0P 
H  r  TJ  P  \’ 

i  E  VC 


B-3 


R-2298 


J 

0 


9 


11 


) 


U  ~C 

O 

o 

•  cc 


function  etc?) 

COMMON  S,M,GGGG 
CG.MVON/COMFT/Dv.(  3  ),SC,  AMD 
CQMyO\’/CJs'l/I.Cl,LC2,LC3,LC*t,LC5 
LC2=LC2  +  1 
FV=V. 

DEL=7*T/(2.*3G) 

Z=DEL*1. 

Z1  =  Z*1 . 

7.C  =  Z*Z 

DEL2  =DEL"DEL 
IF(DZL.GE.3.39)  SO  TQ  11 
I F ( DEL. LT. 1  *  1-7  )  GO  TQ  9 
Q  1  =  QC.  2  (  DEL,  1 ) 

G2  =  G'i2  (  DEL,  2  ) 

G2=G22(2EL,3) 

aFRl=QC(2)/(2.0EL)*(G2-(l.>i/EL  )*C1) 
uFR2=DG(3  )  /  (  2.«-QZL)*(Q3-(  l.+CEL)*G2) 
FT=1./(5G*(2.+DEL)  )*(  (2.*-DEL)/2.*Ql  + 

1  2.  *D:j(2)*((-2.-t-DEL)*Q2  +  (2.  +  DEL-DEL2*(3.+-DEL))*rQl)/ 

1  ((2.f DEL) *DEL)*GPa2-(1.+DEL)*GPR1 ) 

RETURN 

CONTINUE 

FT=l . / ( 2 . *SQ)*((Ey*EM-.75)*ALQG(DEL)+2./DEL) 

RETURN 
CONTINUE 
Q 1 =QG l ( Z, 1  ) 

Q 2=  QQ 1 ( Z, 2  ) 

Q3  =  QG  1 (2,3) 

wPl=DQ(2)/Zl*CQ2-Z*Ql  ) 

GP2=DQ(3)/Z1*('J  3-2*32) 

FT=1./(SU*Z1 )*(Z1/2.*Q1+2.*DC(2)*(((Z-3.)*Q2- 
1  (ZQ*Z-4.*Z*1.  )*Ql)/(Zi-l.)+QP2-Z»aPl)) 

RETURN- 

END 


FUNCTION  ?  T 1 ( ? ) 

COMMON  S , M ,  G 

C0MM0N/C3M1/LC1,LC2,LC3,LC4,LC5 

CCMMCN/COMXX/XX 

LC1=LC1*1 

FT1=(XX-T-G)*FT(I) 

RETURN 

END 


FUNCTION  F T 2  ( T  ) 

COMMON  S , M , G 

CCM.VCN/C3‘!l/LC1,LC2,LC3,LC4,L.C5 

COM  ?•'.  0  .V  /  C  U  .’•!  X  X  /  X  X 

LC1-LC1+1 

FT2=(XX«-7-G)*FT( D 

RETURN 

END 
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PROGRAM  3 

Computes  weighted  average  at  square  of  calculated  axial 
veloci ties. 


FORTRAN 

VARIABLE 


S 

M 

X 

LI 

L2 

Input  same 
as  4  A 

B 

AERR 

RERR 


ANALYTIC 

VARIABLE 


DESCRIPTION 


s 

m 

x  Location  of  propeller  plane 

Determine  the  two  values  of  x  to  be 

[  ° 
used  (both  values  are  same) 

Lower  limit  of  outer  integral 

Upper  limit  of  outer  integral 

Absolute  error  tolerance  for  oyter  i ntegra 1 

Relative  error  tolerance  for  outer  integral 


Output  Cl 


Denominator  of  Equation  ( 1 8 ) 
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FUNCTIJN  FXU) 

F.ATInSAL  F 

c  c ; ;  v  o  n  s  /  v  /  o  -.j  *j  v> 

ccMycN/CG*,';i/cj(5«?'j/3)/ aq(3)/ A^Q/SaO),  GG  ( 3  ), 
Aia(3)/A3a(3),Ea(3)/Ela(3)/33a(3)/FAC(3) 

CCMV 3?l/ CQ:i S 2 /  AO  22  (  6 , 3  ) /  BQ 22  <  6 , 3  ) 

CC‘!uO*!/CDMF  T/DO  (  3  )  /  S-.  ,  AfV 
COMMON/C3*'.F/RRRR/CL/X/X0 
CCMMOS/CDMF X/XX0(4),L1/L2 
CCyy.QN/COMl/LCl/LC2/LC3/LC4/LC5 
IF(ASS(P).CT.l.E-d)  GO  TU  11 
FX=t. 

GO  TO  22 
CONTINUE 
PPRR=R 

PI=3. 1415927 
A  =  -l. 
t  =1 « 

AERP=J. 

RERP=. 30031 
XP=XXD(L1) 

CGI =DCACRE (F,A,3, AERR/RFRR^ERR/  IER) 
SKl=CGl*(l./(d.*r I*3-SQPT(R*3))  ) 

X0=XX2(L2) 

CG2=DCADRE(F,A,?, A ERR/ KERR, ERR, TER) 
SF2=CG2*(l./(3.*?I*S*SwRT(R*S)>) 

F X=R*SK1*SK2 

CONTINUE 

RETURN 

END 

FUNCTION  F ( G ) 

COMMON  S/M,GGGG 

COMMG!,7CCMFT/DQ(3),SG,  AMZ 

CCMyC'I/COMF/R/CL/X/XJ 

CCMM0N/C0M1/LC1,LC2,LC3,LC4,LC5 

LC1=LC1*1 

C*l. 

SK=(X-G)*( X-G) 

RF1=(R-S)*(R-S) 

KF2= ( R+S ) * ( R+S ) 

BKP=4*R*S*DQ(2 )*{ S*S-R* R-3X  )  / ( ( 3K+RK1  )* ( 3K  +  RK2  )  ) 
Z=(3K>R*R+S*S)/C 2*R*3) 

H  =  l. 

IF  ((-C.LT.G). AND.(G.LE.XJ))  PM=H 
IF  ((X2.LT.G). AND. (G.LT.C  ) )  P M= H*(C-G ) / ( C-X J ) 

IF  (Z.LT.1.09)  GO  TO  51 
QlsQQKZ/1  ) 

32*001 (Z, 2) 

GO  TO  m 
CONTINUE 
DEL=Z-1 . 

Q1=CQ2(DEL, 1 ) 

Q2=QQ2 ( DEL/ 2 ) 

F  =  PM*(-Q1*3KP*(Q2-Z*'*1  ) ) 

RETURN 

END 
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PROGRAM  4 

Computer  h  using  least  square  fit 
axial  velocities. 

FORTRAN  ANALYTIC 

VARIABLE  VARIABLE 


between  measured  and  calculated 


DESCRIPTION 


s  s 

M  m 


LI 

Input  L2 

same  as  3  ^ 

3 

AERR 

RERR 

AAAA 


Location  of  propeller  plane 

Determine  the  two  values  of  x  (to  be 

[  ° 

used  (both  values  are  the  same 

Lower  limit  of  outer  integral 

Upper  limit  of  outer  integral 

Absolute  error  tolerance  for  outer  integra. 

Relative  error  tolerance  for  outer  integral 

Denominator  of  equation  (18)  (from  program  3) 


XXO(Ll) 

[Two  values  of  x  used 

XXO ( L2 ) 

o 

Cl 

Numerator  of  Equation  (18) 

AAAA 

x/C 

HHHH 

h 

h  to  give  least  square  fit 

CCLL 

CL 

computer  from  H  obtained  by  this 
program  from  Equation  (13) 

6-13 
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Z  X 


L 


S  ,  ^ / G  j jj 


CC  '.'  j.'i/Cw  *  il/ ^  .  ( 
1  A 1 C  (  3  ) ,  A  3  1  (  ?  ) ,  E 

CCMVj,.’/C3',/l2/  A  ..2 


5?J,3),a:(3)/A22/?3(3),G:;<3) 
U(3),31C(3),S3C(3),FAC(3) 
2(6,3 ),BL22(6,3) 


/ 


CC:'VC?;/CC'<?T/DL  (  3  ),SL, 
CCMMuN/Cj::?  /RRRR,CL,X,X3 
CCMVQN/C3'-'FX/XX.3  (  4  ),  LI ,  L2 
CCMVav/CU^^T/ROV  (9  )  /  ULV.  (9  ) 
CC:!YJN/CCM1/LC1,LC2,LC2,LC4,LC5 


LATA  X  XT /-!.,-. 5,. 7., .6/ 


CAT A  RF//. 1111111, .2222222, .3333333, . 4444444/ 
1  . 5555555,. 6 666666/. 77777 77/. 33 38383/1./ 

LATA  UUY/.  27  31,  .1548,.  2398,.  3*337,.  3316, 

1  . 3496,. 3532, .  32s6, .2975/ 

CALL  ULRSSTCl/LEVL) 


SEAC(5,9J1)  3  ,  M , X , L 1 , L  2 
?CRWAT(F, I,r,2I  ) 

F.  EAC(5,9J2)  A/S/Ac. i  R  ,  P  E  R  R 
FORMAT (4F) 


REA2(5,S«31)  A  A A A 
?I=3. 1415927 


DC  111  1=1,9 

r  r  ( i ) = s  *  r  r  y.  ( i ) 

CONTINUE 


SQ=S*S 


CALL  QCON'ST 

C1=FCADRE(FX,A,B, AERR, SERR, ERR, IER ) 
IF ( AH3( AAAA) .LT. l.E-3)  A  A  A  A  =  1 . 
KHHH=C1 /AAA A 

CCLL=.25*HHHii*(3.+XXid(Ll)) 
WRITE(o,933)  XXtKLl  ),XX0(L2),C1 
W  P 1 7  E  (  6 , 9  0  3  )  A  A  A  A, HHHH,CCLL 
F0R*AT(1X,3E14.6) 
e::d 
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J 

") 


1 


1  1 


51 

1ZU 


FUNCTION  EX  (2) 

£  x  r  e  r  \  4  l  ~ 

w c : . v. o n  ^  / vv uujo 

cnv.»:j‘i/cay  U  /  C  .  (  5  2  2  ,  3  )  ,  AO  (  3  )  ,  A  2  Q,  2  ,1  (  3  )  ,  12  (  3  )  , 
A10(3)/A3^(3)/Z,(3)/Bl'2(3)ybJj(3)/rAC(J) 
CC/.yjN/CGXG:/ A022(C,3  ),CC22(6,3  ) 

C  C  M  M  0  *7  C  C y  F  T  /  3  0  (  3  )  /  S  G  /  A  M  « 

c  a.v.v  3  n  /  coxf  /  rr  r  p  ,  cl  ,  x,  x  z 

C  C  X  V  0  N  /  C  0  y  F  X  /  X  X  2  (  4  )  /  L 1  /  L  2 
CCyyGX/COy 1 /LC1 / LC2 / LC3 / LC4, LC5 
IF(A:i3(R).GT.l.S-S)  GC  TO  11 
FX=J. 

30  TO  22 
CONTINUE 
R  RR  R  =  P 

F  1  =  3. 14159^7 
A  =  - 1 . 

"si. 

AE  RR  =  .' . 
rERR=.UC2Jl 
X  2  =  X  X  J  (  L 1  ) 

CGl  =  DCADfiE(r,A/a/AEPr. /KERR/E  RR/IER) 
SKl=CGl*(l./(d.*:jI*S*S3RT(R*S))) 

5  K  2  =  r  I N  TP  (  R  R  R R ) 

r X=R*SK1*SF2 

CONTINUE 

RETURN 

END 

FUNCTION  F ( G ) 
common  s,y,GGGG 
COMMON/COMc  T/DQ(3  ),SQ,  AM0 
COyyON/COMF/R/CL/X/XJ 
COMM.GN/COyi/LCl  ,  LC2 ,  LC3 ,  LC4,  LC5 
LC1 =LC1 ♦! 

C  =  l. 

3X=(X-3)*(X-G) 

R‘*l  =  (R-S)*(R-5) 

RK2=(R*S)*(R*3  ) 

bK?  =  4*R*S*DQ(2)*(3*S-R*R-8K)/C(  9K  +  RK1  )*(bK«-RK2)  ) 

Z=(9K+R*R*3*S)/(2*R*S) 

ii=I. 

IF  ( (-C.LT.G). AND. (G.LE.X2)  )  PM=H 

IF  ( (XC.LT.G).AND. (G.LT.C) )  PX=H* ( C-G  )/ ( C-X3 ) 

IF  (Z.LT.1.39)  GO  TO  51 
01=031(2/1 ) 

.  2  =  G  .  1  (  Z  /  2  ) 

GC  TO  122 
CONTINUE 
DEL=Z-1 . 

ClsGiZC^tL/l) 

J2=  ,  .2 (  Or  L/  2  ) 
t =P«*( -G1*?NP*( 32-2*31  )  ) 

RETURN 

END 

FUNCTION  FINTP(R) 
ccyyoN/ccyDT/RRM (9 )/Uuy (9  ) 

F  I  \  Tb  =  4  .  *U'JV  (  9  ) 

20  1  1=1/9 
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PROGRAM  5 

Calculates  associated  Legendre  Function  Qm_|/^Z)  for  given 
and  Z.  (Used  as  subroutine  for  all  foregoing  programs.) 
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CC 


J 

1 


cc 

cc 


SUBROUTINE  TC  COS-PJT-  CONSTANTS  DEPENDING  ON  .’•!  A 
INPUT  TO  3  FUNCTIONS 
SUBROUTINE  GC0N5T 

ccxvgn  s,3,gggg 

CGf-.yCN/ CO*'  31  / CG  (  5  0  J,  3  )  ,  AG  (  3  )  ,  A2 C ,  BQ  (  3  ) ,  G3  (  3  )  , 
AlQ(3),A3G(3),EG(3),BlQ(3),o3C(3),FAC(3) 

C  O.MXQ ‘.’/COM  02/A.122(6, 3), EG22(o, 3  ) 

CO.MXGN/COMFT/L)u(3  ),SG,  A*0 

CONSTANTS  FOP  1ST  Q  FUNCTION 


) 


) 


2 

1 

CC 

:c 

cc 


J 


O 

o 


1 1 
21 


) 


40 

3 

CC 
>  CC 


3 


A2G=1. 7724539 

f  v=  x 

DC ( 1 ) =FM- . 5 
DG(2)=F  5 

D  G  (  3  )  =  F  M  ♦  1 .5 
00  1  1=1,3 
AG(I)=(DG( I )  ♦  2 . )  /  2 . 

3G(l)=(DG(I)*l.)/2. 

GO ( I ) =0Q ( I  )  +  1  .  5 
A 1G ( I )  =  DO ( I )  +  l. 

A3*  I)=D J( I ) +1 . 5 
EC(1)=C-;(I  >*1. 

CALL  G.yy'!A(AlG(I  ),31G(I  ),IER1) 

CALL  G  M  M  A  (  A  3Q  ( I )  ,  B  3G  ( I )  ,  I  £  R  3 ) 

r AC(  I ) =B IQ ( I )*A2G/( (2.*»SG( I ) )*B3G(I )  ) 

CG(1,I)=A3(I)/GG(I) 

DO  2  L=2,5D0 

CC(L,  I )=C3(L-1, I )*(< AG( I J  +  L-l . )/(GQ( I  )  +  L-l.  )) 

CONTINUE 

CONTINUE 

CONSTANTS  FOP  2ND  0  FUNCTION- 

DO  3  L  =  1 ,  3 
\22  =  QG( L)+  .5 
CG22=0. 

IF(N22.E*.D)  GO  TO  21 
DO  11  1=1, N  2  2 
CG22=CG22*1 ./(  2.  *1-1.  ) 

CONTINUE 

CONTINUE 

AG22(l,L)=S./2.*ALOG(2. )-2.*CG22 
5  G  2  2 ( 1 ,  L ) =~ . 5 
CO  40  1=1,5 
NP  =  I  - 1 

BG22(  I  +  1,L)  =  L<  322(1, L)*(N  22**2  •l./4.-N?*CN’P  +  l))/ 
1  (2.*(N?*1)**2) 

AG22 (  dl,L)  =  A022( I,L  )*(S22**2-1 ./4.-NP* (NP  +  1) )/ 
1  (2.  *(  NP*1 )  **2)-BQ22(  I,L)*(2.  *(  N'22**2-l./4.  )♦ 

1  (  N  P  ♦  1 )  )/(2.*(NF*l)**3) 

CONTINUE 

CONTINUE 

CALCULATION  OF  ,US£D  PY  FUNCTION  FT 

A  y  l  =  o . 

DO  3  X=l,v 

ami  =  ami«-i./(2.*k-i.  ) 

CONTINUE 

AVU=2.5*ALOG(2.)-2.*A::i 


RETURN 

end 
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first  ;  function 

^  U  N  C  ?  I  _  N  i.lC/J 

cc:-.u j’-’/cg':',  l/c  j ),  a„  c  ),  ^  j(  3 ),  c  j(  j  ), 

IK  (  3),  A  9  ;(  3),  £  ;■(  3  ),  51  0(3  ),  c  Jv.  (  3  ),  EAC(  3) 
comeon/covet/d  ,(3 ) ,sg,  am.- 
C C -vi y  2  N/ CJ'U  /  LC1  ,  LC2,  LC3  ,  LC4 ,  L  C5 
LC3  =  LC3+  1 
ZF=l.y (Z*2) 

F  =F  AC (L)/ (Z**EO(L) ) 

DCJ=9G(L)*F*ZF 
DC  2  3  1=2/  5 iV 
F  L  I  =  I 

T  =  C  2  (  I-l/L)*Dr>D 
F  =  r  +7 

IF(  iPS(T)  )  GO  TO  5 

CDD=D9D*ZF*  (  (3G  (D  +  FLI-l  .  )/FLI  ) 

CD'.’T  I  VUE 
CONTINUE 
Q  w  1  =  F 
RETURN 
EMC 

SECOND  ,/  FUNCTION 
FUNCTION  QU?(DEL,L) 

COM  V0V/C3''  22/ A  .22(0, 3),  SO  22  (0/3  ) 

com  von/ com  1/  lci ,002,003,0*04,005 

1.09  =  0044-1 
SUM  1 =3 . 

3  U.*!  2  =  ^ . 

T"STJ1=«5. 

TESTS 2= 8. 

DC  40  1=1,5 
N  P  =  I  - 1 

Sl’Vl  =SUV1*AQ22(  I ,  L  )  *  DEL  **.NF 
3UM.2  =  5U/2+5U2  2  (  I,L)  *  CEL**  UP 
TEST1=A9S( 3UM1  ) 

T£ST2=APS( SUM2 ) 

IF(AHS(TEST1-TE3TC1).GT.0. J338J31  )  GO  TO  40 

IF( AES(TE3?2-TEST02) .LT. 0.3000301)  GO  TO  50 

TEST31=7ES71 

TES  T J  2=  TES  T2 

CONTINUE 

CONTINUE 

QQ2=SUHU ALOGC DEL) *3UM2 

RETURN- 

END 
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Figure  1.  Radial  Variation  of  Shaft-Frequency  Wake  Harmonic 
for  NSRDC  Model  4602  (from  Hadler  &  Cheng  SNAME 
Trans.  ,  Vo  1 .  73 ,  '965) 
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